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ABSTRACT 

Q I We investigate the 21 cm emission from the epoch between reionization z^. 

CN I and the Gunn-Peterson transparency zqp- According to the lognormal model of 

Q I the thermal history around reionization, hydrogen clouds m. z^ > z > zqp are 

O \ hot and a predominant part of baryonic gas is ionized, but still opaque to Lya 

CN ' photons. Therefore, 21 cm emission is a distinctive characteristic of this epoch. 

I We show that the 21 cm emission comes from both uncollapsed and collapsing 

CS I hydrogen clouds. The spatial distribution of the brightness temperature excess 

t;j- ' 8T\, is highly non-Gaussian. It consists of spikes with high bT^, and a low ^T^ 

J^ ■ area between the spikes. The field has the following statistical features: (1) the 

\^ • one-point distributions of 6Tb are described approximately by power-law tailed 

^ ' probability distribution functions; (2) the nth-order moment of 5Th is increasing 

O ■ much faster with n than that of a Gaussian field, but slower than that of a 

1-^ ■ lognormal field; (3) the scale-scale correlation of 6Tb field is significant for all 

I ' scales larger than the Jeans length of the gas. These features would be useful 

Vh ' for distinguishing the 21 cm emission of the early clustering from the noise of 

^ ■ foreground contamination. 

C^ ■ 



% 



Subject headings: cosmology: theory - intergalactic medium - large-scale struc- 
ture of the universe 



1. Introduction 

The universe becomes transparent to Lya photons (the Gunn-Peterson transparency) 
at redshift zgp — 6.5 (Fan et al. 2002). Yet the polarization of the cosmic microwave 
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background (CMB) indicates that the Compton scattering optical depth to the CMB is as 
high as Tg ~ 0.17 ± 0.04 and that the reionization of the universe probably occurred at 
11 < Zr < SO (Kogut et al. 2003). Most theoretical attempts to discriminate between the 
redshifts Zcp and Zr concluded that the history from first-generation star formation to the 
Gunn-Peterson transparency lasts a long time (Cen 2003; Haiman & Holder 2003; Whyithe 
& Loeb 2003; Holder et al. 2003; Onken & Miralda-Escude 2003). 

This problem has recently been studied by using the lognormal (LN) formalism (Bi et al. 
2003; Liu et al. 2004). Considering that the Jeans length changed greatly before and after 
the reionization, a self-consistent calculation with the LN model yields a general picture that 
the redshift Zr should be significantly larger than zgp, where Zr means the era that ionized 
regions are fully overlapped. The thermal history of the universe around the reionization 
can be roughly divided into three epochs: (1) the cold dark age z > Zr, (2) the hot dark age 
Zr > z > Zgp, and (3) the bright age z < zqp- With the cosmological parameters given by 
WMAP and COBE, the hot dark age lasts about Zr — zqp — 10. 



In the epoch Zr > z > zqp, the fraction of ionized hydrogen fnjj = nHii/ijiHi + nnii), 
where uhi and nnii are the number densities of neutral and ionized hydrogen atoms, is 
already significant, while the mass fraction of neutral hydrogen fni = nHi/ijiHi + uhii) is 
still much greater than 10~^, so that the universe is opaque to the Lya photons. Moreover, 
the temperature of hydrogen gas in the period Zr > z > zcp is higher than that of the CMB. 
Therefore, in this epoch, the 21 cm emission given by the hyperfine structure of the ground 
state l^Si/2 of neutral hydrogen would be significant. This motivates us to investigate in this 
paper the statistical properties of the redshifted 21 cm emission from the epoch Zr > z > zgp- 

Many calculations on the redshifted 21 cm radiation from the early universe have been 
done based on various models of the thermal history around reionization (Ciardi & Madau 
2003; Furlanetto & Loeb 2002, 2004; Iliev et al. 2002, 2003; Tozzi et al. 2000; Zaldarriaga et 
al. 2004; Furlanetto et al. 2004). These works have generally focused on the power spectrum 
of the fluctuations of the 21 cm emission. In this paper, besides the power spectrum, we 
address the higher order statistics and non-Gaussianity of the spatial distribution of the 21 
cm emission. Since the 21 cm emission at high redshifts is relevant to the earliest gravitational 
clustering and star formation of the universe, the field of the 21 cm emission must be non- 
Gaussian. It would probably be the earliest non-Gaussian field if the initial perturbations of 
cosmic mass density field are Gaussian. 

The non-Gaussianity of 21 cm emission field is not only theoretically important but 
also observationally useful. Radio observation would be able to detect the cosmic 21 cm 
emissions if its fiux is of the order of ^ 0.1 mJy or higher and angular scale is of the order 
of arcminutes (Morales and Hewitt 2004; Pen et al. 2004). However, it is a challenge to 
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identify the redshifted 21 cm emission from the recombination epoch, since the noise from 
the artificial radio interference in the VHF band is serious. The non-Gaussian statistical 
features would be helpful to draw the information from the noisy observations. 

This paper is outlined as follows. §2 briefly describes the basic results of the reionization 
and thermal history of hydrogen clouds in the LN model. §3 presents the simulations of the 
baryon field and the 21 cm emission from neutral hydrogen in the epoch of z^ to zqp. 
The basic properties of the 21 cm emission from diffused intergalactic medium (IGM) and 
collapsed halos are also addressed. §4 presents a statistical analysis of the bT^ fields, including 
the first-, second- and high-order statistics. Finally §5 gives the discussion and conclusion. 



2. Evolution of hydrogen clouds around reionization 

2.1. Hydrogen gas distribution in LN model 

The LN model assumes that the mass density distribution p(x) of hydrogen gas at 
redshift z is given by an exponential mapping of the linear distribution of the dark matter 
mass density contrast 5o(x), which is smoothed on the Jeans length scale of the gas at the 
redshift z (Bi 1993), i.e., 

p(x) = poexp[5o(x) - (7o/2], (1) 

where po is the mean density and o"o = (^o)^ i^ ^^e variance of the Gaussian mass field 5o(x) 
on the scale of the Jeans length of the gas. The probability distribution function (PDF) of 
the density field p(x) is, then, lognormal: 



p{p/ P) = 1 1= exp 

(p/p)o'oV27r 



1 / ln(p/p) + aoV2 

2 V ^0 



p > 0. (2) 



The dynamical arguments of the exponential mapping and LN PDF of both the dark matter 
and hydrogen clouds have been studied by, e.g.. Coles & Jones (1991), Jones (1999), and 
Szapudi & Kaiser (2003). 

The first characteristic of the LN model is that the clustering of hydrogen clouds is 
mainly dependent on the variance o"o of the linear fluctuations on the scale of the Jeans 
length, while the clustering of dark matter basically is independent of the properties of 
gas. Therefore, the LN model is effective in describing the effect of the thermal status of 
hydrogen gas on the density and velocity field of the baryonic gas. This has been shown 
in the successful application of modelling IGM Lya forests of QSO's absorption spectra in 
redshift range 2 < z < 5 (Bi 1993; Bi & Davidsen 1997; Feng & Fang 2000). 
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The second characteristic is that, compared to a Gaussian or a normal distribution, the 
PDF of an LN model described by equation (2) has a prolonged tail at the high-density end 
and that the probability of high- density events p/p ^ 1 is very sensitive to o"o- This property 
has already been found to be useful in explaining the intermittence of the transmitted flux 
of QSO's Lya absorption spectrum (Jamkhedkar et al. 2000, 2003; Pando et al. 2002). The 
long tail can directly be seen with (1) the cumulative mass fraction M(> p/p)., which is the 
fraction of mass in regions having mass density larger than a given p, and (2) the volume 
filling factor V{> p/p), which is the fraction of volume with density larger than a given p. 
They are 

M[> (pIp)] = f M':)dx = ierfc (i^J^ - ^) , (3) 

Therefore, even when o"o — 1, the mass fraction of high-density events is already significant 
and their volume filling factor is very small. 



2.2. The clustering of hydrogen clouds at redshifts Zr and zqp 

In the context of the LN model, the evolution of hydrogen clouds is mainly dependent on 
the Jeans length of hydrogen gas. If hydrogen gas is uniformly distributed, its Jeans length 
is given by \j = fs(7r/G'p)^/^, where Vg is the sound speed of the gas. For cosmological study, 
it is convenient to use a comoving scale xj = \j/27r = {1/ Ho)[2'jkBTm/'ifxmpil{l + z)]^^"^, 
where Tm is the mean temperature, p the molecular weight of the gas, Q is the cosmological 
density parameter of total mass, and 7 is the ratio of specific heat. 

Primordial baryons, created at the time of nucleosynthesis, recombine with electrons 
to become neutral gas at 2; ~ 1000. Before z ~ 200, the residual ionization of the cosmic 
gas keeps its temperature locked to the CMB temperature (Peebles 1993). After z ~ 200, 
the gas cools down adiabatically because of the expansion of the universe. Assuming the 
adiabatic index 7 = 5/3, i.e., hydrogen temperature T oc p^^^, the comoving Jeans length 
evolves as ~ 0.1 x (1 + zY^'^ h^^ kpc. During reionization, the gas is heated by UV ionizing 
photons from a low temperature to ~ 1.3 x lO^K. This leads to the increase of xj at Zr. 

Assuming Zr — 18 (Kogut et al. 2003; Liu et al. 2004), the evolution of the Jeans 
length is sketched in Figure 1, which gives the relation between xj and the cosmic scale 
factor 1/(1 + z). The sudden increase of xj at Zr shown in Figure 1 is certainly unrealistic, 
which comes from the assumption that the gas temperature increases instantly by a factor of 
about 10^ at Zr- More realistically, the sharp changes of xj at Zr and zcp should be replaced 
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by a softened transition. Nevertheless, the physical status of the gas before and after Zr is 
properly sketched in Figure 1. 

Corresponding to the evolution of xj at Zr-, the variance ctq, which is the rms mass 
fluctuation on the Jeans length [^pj p)xj oc k^P{}i)\k=\ixj-, also undergoes a zigzag evolution 
around Zr- Figure 2 shows o"o versus 1/(1 + z) calculated with the xj of Figure 1. To 
calculate (Tq, the power spectrum P{k) of the linear mass density perturbations is taken to 
be the spectrum of the low-density fiat cold dark matter model (LCDM), which is specified 
by the density parameter Qq = 0.3, the cosmological constant ^a = 0.7, and the Hubble 
constant h = 0.7. The linear power spectrum P{k) is given by the fitting formula given 
by Eisenstein & Hu (1999). The increase of o"o with 1/(1 + z) shown in Figure 2 is due to 
the gravitational linear growth. The variance at small redshifts {z < 7) shown in Figure 2 
actually is the same as that given by Bi & Davidsen (1997). The zigzagged feature at z ~ 3.3 
of Figure 2 is due to Hell reionization (Theuns et al. 2002). 

An important feature of o"o shown in Figure 2 is that the curve of (Tq from A' to B' 
almost repeats that from A to B. Thus, the evolution of the fraction of mass [eq.(3)] in the 
period from A' to B' is similar to that from A to B. The period from A to B corresponds 
to an evolution of weak-to-strong clustering of the baryonic gas, and the period from A' to 
B' also corresponds to an evolution of weak-to-strong clustering of baryonic gas. Therefore, 
the formation rate of collapsed hydrogen clouds at A' is significantly lower than that at B. 
This leads to a suppression of clustering and star formation just after Zr- 

The dynamical picture of the suppression of clustering can be explained by the following 
negative feedback mechanism. During the reionization Zr, hydrogen clouds are heated by the 
ionizing photons so that the mean temperature of hydrogen clouds increases by a few orders. 
All the irregularities originally formed in cold gas on scales smaller than the new Jeans length 
would be smoothed out by the heated gas, and the variance of the fluctuations of hydrogen 
cloud distribution drops back below unity. Although dark matter halos continuously collapse 
before and after Zr, hydrogen clouds on scales less than the new Jeans length will stop 
collapsing. As a result, star formation will be slowed down or may be even halted in these 
halos after Zr- The level of the star formation suppression is adjusted by the following 
mechanism. Once the star formation rate declines, the UV background produced by the star 
formation also becomes lower, and so do the temperature and entropy states of hydrogen 
gas, which make it easy for further star formation. On the contrary, if the star formation 
rate increases, the UV background produced by the star formation also increases, and then 
more hydrogen gas will be heated to higher temperature and higher entropy states. This 
finally yields more suppression in the rate of star formation. 

When the variance (Tq grows again to unity around B' because of the increase of the 
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potential wells of dark matter, there is once again enough to have collapsed hydrogen clouds, 
star formation, and ionizing photons. This finally gives rise to the Gunn-Peterson trans- 
parency. Therefore, the long-lasting period from z^ to zqp is due to the zigzag of gq around 
reionization z^. The difference z^ — zqp is found to be equal to ~ 10. It is weakly dependent 
on, or stable with respect to, the following parameters: (1) the threshold of the collapsed 
hydrogen clouds pjp and (2) the mean ionizing photons released by each baryon in the col- 
lapsed objects (Liu et al. 2004). This stability is due to the above-mentioned adjustment 
mechanism. In this sense, the redshift of reionization z^ ^ 18 actually is determined by the 
redshift of the Gunn-Peterson transparency zgp (Liu et al. 2004). This result is also in good 
agreement with the observed optical depth t^. to the CMB photons. 



2.3. Temperature of hydrogen clouds in the epoch of z^. to z.^ 
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According to the above-discussed scenario, in the epoch z > Zr, hydrogen gas is cold and 
neutral. It may contain some individual ionized spheres around the first-generation stars. 
Hydrogen clouds generally are opaque to Lya photons. When z approaches Zr, the increasing 
individual ionized regions overlap. At Zr, the individual ionized regions fully overlap in the 
sense that the universe becomes transparent to soft X-ray photons. Thus, after Zr, hydrogen 
clouds are heated by soft X-ray photons and the star formation is suppressed. The fraction of 
ionized hydrogen {HI I) will no longer increase, or will even decrease after z^-. Consequently, 
hydrogen clouds are still opaque to Lya photons. 

It has been shown that if the photoionization heating is the major heating process of 
the hydrogen clouds, the mean gas temperature is always of the order of ~ 10^ K, weakly 
dependent on the intensity of heating photons (Black 1981; Ostriker & Ikeuchi 1983). This 
result is further supported by hydrodynamic simulations of baryonic gas in the early universe 
(e.g., Hui & Haiman 2003). They showed that the mean gas temperature asymptotically 
approaches 10'^ — 10^ K. This temperature is not determined by the intensity of the heating 
photons, but is completely given by the shape (or hardness) of the spectrum of the heating 
photons. The mean temperature is higher for harder photon spectra, such as soft X-rays. 
Thus, although star formation is slowing down after Zr-, the mean temperature will still be in 
the range 10^ — 10^ K. Therefore, it is reasonable to assume that in the epoch of Zr to zgp-, 
the mean gas temperature can be described as T^ — Tq[{1 + z)/{l + Zr)]^ , where Tq ~ lO^K 
and the index tj is of the order of 1. In Figure 3, we plot the mean temperature evolution 
for the two cases of (Tq,?]) = (0.8 x lO^iT, 2) and (1.3 x IQ'^K,2). In Figures 1 and 2 we 
also show, respectively, the temperature effect on the Jeans length xj and variance ctq in the 
epoch Zr > z > Zgp- 
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If hydrogen gas is polytropic, gas temperature at the point with mass density p is given 
by T = Tmip/pY'^^- In this way, we can calculate the thermal status of hydrogen from 
the density distribution. In the LN model, the relation T — p and the index 7 actually are 
given by fitting the hydrodynamical simulation (Bi & Davidsen 1997). In this case, shock 
heating is partially considered. However, if gravitational shocks are strong, the temperature 
of hydrogen gas is basically multiphased, and one cannot describe the T — p relation by a 
single equation. To estimate the effect of strong shock, we have to replace the temperature- 
density relation T = Tm{p/p)^^^ with P{T;p), the PDF of temperature at a given p. Such 
a PDF has been recently derived by using a hydrodynamic simulation code that is very 
effective in capturing shocks (He et al. 2004). He et al. (2004) found that strong shock 
events are very rare at 2; > 4, i.e., that the events leading to the deviation from the relation 
T = Tm{p/p)"'~^ are negligible at 2; > 4. 



3. Samples from the LN simulation 

3.1. Simulation of the hydrogen gas distribution 

We produce simulation samples of spatial distribution of hydrogen gas p(x) by the LN 
model. In order to quickly grasp the features of these distributions, we only concentrate on 
samples of one-dimensional distribution. The details of the simulation procedures have been 
given in Bi & Davidsen (1997) and Bi et al. (2003). A brief description is as follows. 

We first generate the one-dimensional density and velocity distributions in Fourier space, 
So{k, z) and v{k, z), which are two Gaussian random fields. Both So{k, z) and w(/c, z) are given 
by the power spectrum, Po{k) as follows (Bi 1993; Bi et al. 1995) 

6o{k,z) = D{z){u{k)+w{k)), (5) 

TT 

v(k, z) = F(z)—ika(k)w(k), (6) 

c 

where D{z) and F{z) are the linear growth factors for the So{x) and v{x) fields at redshift 

z. The w{k) and u{k) fields are Gaussian with the power spectra given by 

POD 

PUk) = a-' / Poiq)2iTq-'dq, (7) 

poo 

Pu{k)= Po{q)27iqdq-P^{k), (8) 

Jk 



where Po{k) is the power spectrum of the three-dimensional field 5o(x). The functions a{k) 
in equation (7) is defined by 

The power spectrum Po{k) is taken to be (Bi & Davidsen 1997; Bi et al. 2003) 

Obviously, Po{k) is also a function of 2; via the redshift dependence of xj{z). Thus, for 
a given Pouik) and xj, one can produce the distributions 6o{k,z) at the grid points ki, 
i = 1,2,.. .,N in one-dimensional Fourier space. The spatial distributions 6o{x,z) can be 
obtained by using the fast Fourier transform. Since the velocity follows the linear evolution 
longer than the density, we can use the linear v{k, z) and its Fourier counterpart as good 
approximations of the velocity field at high redshifts. In Figure 4 we plot a typical realization 
of the gas mass density fields at redshifts 7, 10, and 15. The simulation size shown in Figure 
4 is 50 h~^ Mpc in comoving space. The total number of pixels is 2^^ = 16, 384. The pixel 
size is then ~ 3 h~^ kpc, which is less than the smallest scale of xj at z < Zr. Therefore, 
the samples are qualified for studying the non-Gaussian features of the hydrogen clouds on 
small scales. 

As analyzed in §2.3, the effect of strong gravitational shock at high redshifts is negligible 
and the temperature of hydrogen gas can be calculated by the relation T = Tm{p/p)'^~^- We 
then have the temperature field of the hydrogen clouds. With the mass density and temper- 
ature fields, one can find the neutral hydrogen mass density fields with the photoionization 

equilibrium equation as 

a{T) 

^""^ a(T) + r(T) + J/ne' ^ ^ 

where Oi{T) is the recombination rate, T{T) is the collisional ionization rate, Ue is the number 
density of electrons, and J the rate of photoionization of hydrogen (Black 1981). Although 
equation (11) is reasonable for determining fni of the IGM aX z < zgp, there are two points 
should be addressed when applying it to the epoch Zr — zgp- 

First, equation (11) is applicable only if the hydrogen clouds are optically thin for 
photons of J. However, hydrogen clouds oi Zr- > z > zgp are optically thick for Lya 
photons and may also be so for photons of ~ 13.6 eV. Therefore, in our case, the fiux J in 
equation (11) is not for Lya, but for soft X-ray photons. As discussed in §2.3, in the epoch 
Zr > z > Zgp, hydrogen clouds are transparent to soft X-ray photons, which are the major 
sources of photoionization heating. This picture is consistent with the fact that the photon 
spectra from the first-generation (massive) stars are generally hard. 
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Second, in the case oi z < zgp, the parameter J is determined by fitting it with the 
mean transmitted fiux of Lya photons (e.g., Choudhury et al. 2001). Obviously, no such 
observations can be used to constrain the parameter J in the epoch oi Zr > z > zgp- 
However, J/Ug can be adjusted by fitting it to the mean fractions of HI or HII ^ which are 
required by the optical depth Tg. Using this method, we can produce samples of fui field. 
Both the models of (Tq, ?7) = (0.8 x 10'^ -ft', 2) and (1.3 x 1Q*K, 2) are consistent with the optical 
depth Te from the WMAP. 



3.2. Simulation of the 21 cm emission of hydrogen clouds 

With the sample of gas mass density, temperature, and fni distributions, we can calcu- 
late the 21 cm (z/q = 1420 MHz) emission of neutral hydrogen (HI). This emission at redshift 
z is determined by the difference between the spin temperature of neutral hydrogen, T<j(z), 
and the temperature of the cosmic microwave background, Tcmb{z) = 2.73(1 + z) K. There 
are two mechanisms leading to Ts{z) > Tcmb, collision and radiation background, which 
gives (Field 1958, 1959) 

I + Vc + Vl 
where Tc is the temperature of hydrogen gas, Tl is the temperature of Lya photons, and 
yc and yi are the collision and radiative excitation efficiencies, respectively. Since hydrogen 
clouds are optically thick to the Lya photons, it is reasonable to assume that Lya photons 
are in approximate thermal equilibrium with the gas. Thus, one can take T^ to be the gas 
temperature. 

Considering self-absorption, the brightness temperature of the 21 cm radiation at red- 
shift z is determined from the radiation transfer equation as 

Ttiz) = T,„,(z)e-^(^) + r ' T,(/)e-^(^')rfr(^'), (13) 

Jo 

where the first term on the right-hand side is from the CMB and the second term is from 
neutral hydrogen; t{z) is the optical depth of the 21 cm absorption. When T^ is much larger 
than T, = /iz/q/^b = 0.06 K, t{z) is given by (Wild 1952) 

where factor F{z, z') is the normalized line profile. For Doppler broadening, we have 

1 -( jI^I-V 

F(z z') = e vHi+^)J CIS) 
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where b = {2kBT/mc^Y^'^. For gas with temperature < 10^ K, we have b < 10~^, and 
therefore the integral of equation (14) lasts for a very narrow range z ± Az and Az < 
10~*^(1 + z). Thus, the velocity distortion can be ignored in the first stage. 

If nHi{z) and Ts{z) are not strong functions of z, equation (14) gives approximately 
(Field 1959) 

_ ^hc-AiQnHi{z) 
'^^^^ ~ ?>2TiulkBTs{z)H{zy ^ ' 

Therefore the fields nHi{z) and Ts{z) in equation (16) are, respectively, the HI number 
density and spin temperature Ts{z) smoothed by the window function equation (15). Thus, 
equation (13) yields approximately Tb{z) = TcMB{z)e~'^^^^ + Ts{z){l — e""^*^^^), and therefore 
the observed brightness temperature excess at the redshifted frequency u = uq/{1 + z) \s 

6n{u) = [T,{z) - T,^b{z)] ^ ~ l^ ^ . (17) 

We show in Figure 5 a typical realization of the fields of T^, t{z), and (5Tb at z=7 in the 
temperature model (Tq, rj) = (1.3 x lO'^K, 2). Figure 5 shows that the field of the brightness 
temperature excess, 6Ti,, is highly nonuniform. At most places, 6Ti, is actually zero but 
contains many high spikes, or patchy structures. The average of STb is only ~ 2.3mK, while 
the spikes can be as high as 20 - 40 mK. The rms of the brightness temperature fluctuations, 
{6T^y/^, is ~ 7 mK. For {To,r]) = (0.8 x 10^^,2), we have {6T^y/^ = 9.5 mK. It has 
been argued that the brightness temperature excess would be too small to observe when 
considering the foreground noise contamination. However, the nontrivial features of the 
non-Gaussianity might help identify the redshifted 21 cm emission from the contaminated 
data. 

The non-Gaussianity of the 6Ti, field can also be seen with (d{6T^)^^'^ / dp)dp, which 
is the fraction of {5T^)^^'^ given by regions with hydrogen gas density from p to p+dp. 
Figure 6 plots the integrated fraction function {5T^y^'^{> p) = J^ d{ST^y^'^. It shows 
that about one-third of {5T^)^'^ is given by regions with density 2 < p < 6 and that the 
other two-thirds are from regions of p > 6. Because the p > 6 regions correspond to the 
collapsed/collapsing hydrogen clouds (Bi et al. 2003; Liu et al. 2004), Figure 6 shows that 
both diffused and collapsed hydrogen clouds of the mass field have comparable contributions 
to the brightness temperature fiuctuations of the redshifted 21 cm emission. Therefore, the 
non-Gaussianity caused by the quasi-linear and the nonlinear clustering of hydrogen clouds 
should be significant. 
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4. Statistical properties of 21-cm emissions 

4.1. DWT variables of 5Tb field 

In order to effectively describe statistical features of (5T{, field, we decompose the 5Ti, 
field by the discrete wavelet transform (DWT; Fang & Feng 2000; Fang & Thews 1998). 
Because the DWT modes are spatially localized, each dimension can be treated separately, 
and hence it is easy to match the simulation with the physical condition of the sampling. 
The three-dimensional power spectrum can be effectively recovered by DWT from the data 
in an area with "poor" geometry. Therefore, with the DWT, one can effectively compare and 
test the predicted statistical properties from one-dimensional or two-dimensional simulation 
with various observations (e.g., Guo et al. 2004 and references therein). 

In the DWT scheme, there are two sets of bases given by (1) scaling functions (pj^iix) = 
{x\j,l)s and (2) wavelets ipj,i{x) = (x|j, /)^, where j = 0,1,.. and / = 0, 1...2-' — 1. In 
one-dimensional space with size L, the scaling function 4>j,i{x) is localized in physical space 
lL/2^ < X < {I + l)L/2^ while wavelet ipj^i is localized in both physical space lL/2^ < x < 
(/ + l)L/2-'' and Fourier space tt2^/L < \k\ < {3/2)2tt2^ /L. 

The scaling functions are orthonormal with respect to the index / as 

s{j,l\j,l')s = j (t>,M)(t>u'{^)dx = S^i,. (18) 

Scaling function (f)j,i{x) actually is a window function for the spatial range lL/2^ < x < 
{I + 1)L/2K The scaling function coefficient (SFC) of a field 5Tb{x) is defined as e^y = 
s(j, l\5Ti,) = f 6Th{x)(l)j^i{x)dx, and therefore, the mean oi 6Ti,{x) in the spatial range lL/2^ < 
x<{l + l)L/2^ is 

^Tum = ^^^^ (19) 



''^^''^ ,(//|l)' 



where |1) is a uniform field with field strength equal to unity. 

The wavelets ipj,i{x) are orthonormal bases with respect to both indexes j and / 

s{j,l\j',l')s = / '4)j,i{x)i)y^v{x)dx = Sjj>6i^i>. (20) 

The wavelet function coefficient (WFC) of the STh{x) field is defined as 

ej,i= w{j,l\S) = 6Tb{x)iljj^i{x)dx. (21) 

Since the set of the wavelet bases \j,l)w is complete, the STb{x) field can be expressed as 

5T,(x) = 5^5^e,>,v(a^), (22) 
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where each j runs 0, 1, 2... and / runs 0, 1,...2'' — 1. Therefore, the WFCs ij^i can be used 
as the variables of the field 5Tf,. The ij^i are the fluctuations around scales k = 27m /L with 
n = 2^ and at the physical area / with size Ax = L/2K 



4.2. One-point distributions of ST^^i^j^i) and STi,(^j^i-) 

We plot in Figure 7 the one-point distributions of 6T{,^(j^i) for two temperature models 
(To, T], z) = (1.3 X lO^K, 2, 7) and (0.8 x lO^K, 2, 10). The scales j are taken to be 6, 8 and 
10, respectively, corresponding to comoving smoothing sizes 0.78, 0.20, and 0.05 h~^ Mpc or 
angular resolutions 0'.44, O'.ll, and 0'.028 at z=7, or 0'.41, O'.IO, and 0'.025 at z=10. All the 
PDFs in Figure 7 seem to have similar shapes. A peak at 6Th ~ is given by the low mass 
density (0 < p < 2) areas; Figure 6 shows that the areas with mass density < p < 2 do not 
contribute to 6Ti). Except for the peaks, the PDFs of Figure 7 are flat in the range 6Th^(^j^i) < 1 
mK, and approximately of a power law at 6Ti,^(^j^i^ > 1 mK. The long tail obviously is given 
by the spikes in the STh field (Fig. 5). The power law tail shows a slightly smoothing scale 
dependence. 

We plot in Figure 8 the PDFs of the WFCs, in which ij^i is replaced by STbi^j^i-j = 
{2^ /Ly/'^ij^i since the latter has the dimension of temperature. The parameters used in 
Figure 8 are the same as those in Figure 7. The PDFs shown in Figure 8 are also weakly 
dependent on the temperature parameters, but the tails of the PDFs are more significantly 

—a 

dependent on scale j than in Figure 7. If we describe the tail by a power law oc 6Tbfji\, the 
index a is bigger with larger scales. In other words, the tails are longer for smaller scales. 
This leads to the mean power {STb^j n) being lower for higher j (see next subsection on power 
spectrum). 



In either Figures 7 or 8, the mean of ST^^q^i) or ST^j^i-^ is given by the events corresponding 
to the long tail of the PDFs. That is, the maps of 5Tb or 5Tb(^j^i) are very different from the 
Gaussian noise field. 



4.3. Second-order correlations 

The two-mode correlation of the SFCs {5Tb^[j^i)5Tb^(j^ii)) is similar to the ordinary two- 
point correlation function. Figure 9 presents the correlation function {6Tb^(j^i)6Tb^(j^i')) for 
parameters {To,r],z) = (1.3 x 10*^ K, 2, 7) with smoothing scales j = 12, 10 and 8, which 
correspond, respectively, to comoving sizes 0.03, 0.05, and 0.20 h~^ Mpc, or 0'.007, 0'.028 and 
O'.ll. The two-mode correlation functions are typically of power law (5T;, (^ /)5T;, (^y/)) oc r~'^ 
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in the range 0.2 < r < 1 h^^ Mpc with index 7 ~ 0.6. The central part of the correlation 
function, or r < 0.2 h'^ Mpc, is flat because of the Jeans length smoothing. The tail, or 
r > 2 h~^ Mpc, of the correlation function approaches (5T;, (ji)(5T;, (j^^/)) = (^T^ q- ;))^ = (STh)^, 
which is the mean of 6Th. 

The two-mode correlation function of the WFCs {^j,iij,i'), or {6Tb^(^j^i^6Tb^(^j^i/)), generally 
is diagonal with respect to (/, /'), especially if the initial mass perturbation is Gaussian (Feng 
& Fang 2004; Guo et al. 2004). Thus we have in general 

(5T,,(,,,)<5T,,(,,,)) = P,5,5„ (23) 

where Pj is the DWT power spectrum of the one-dimensional fleld, of which the relation 
with the Fourier power spectrum P{k) has been given in Fang & Feng (2000). Figure 10 
gives the power spectrum Pj for parameters (Tq, 77, z) = (1.3 x 10'' K, 2, 7) and (0.8 x 10'' K, 
2, 10). The power spectra of these two cases are similar. However, if using angular scales, 
they are different. For instance, although both spectra are peaked at j = 7, they correspond 
to angular scale 0'.22 at 2; = 7, or 0'.20 at z = 10. 

The power (5Tb /• n) increases from j = 2 (comoving scale 12.5 h Mpc) to j = 7 (0.39 

- — 2 
h~^ Mpc). This is because the clustering is stronger on small scales. The power {6Tl|r■l^) 

gradually decreases with j when j > 7 (comoving scale less than 0.39 h~^ Mpc), this is 

because of the Jeans smoothing of hydrogen gas. 



4.4. High-order moments of 5Ti, field 

High-order moments are sensitive to non-Gaussianity. We use two standard high-order 
moment statistics: (1) {{6Tb-{6Tb)f''y/''/{{6Tb-{6Tb)f) and (2) {{6Tlljy/''/{{6Th)). If 
the 6Tb fleld is Gaussian, the high-order moment of 6Th satisfles the relation 

On the other hand, if the 6Tb fleld is lognormal, we have 

^-^^^§^ = exp[{l/2){n-iy], (25) 

where a is the variance of ln{6Tb). 

We plot in Figure 11 the high-order moments {{6Tb - ((5T6))2")1/«/((5T5 - (^Tfe))^) and 
{{6T^jy/'' / {{6Tb)) for the parameters {To,ri,z) = (1.3 x WK,2,7) and (0.8 x lO^/sT, 2, 10). 
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As a comparison, we also plot the corresponding Gaussian and LN moments in Figures 11a 
and lib, respectively. Figure 11 shows that the STi, field is neither Gaussian nor lognormal. 
The high-order moments are always significantly higher than a Gaussian field for n > 1, 
but always less than an LN field. In both temperature models, the two high-order moments 
are always quickly increasing with n when n < 5, and then slowly increasing when n > 
5. Therefore, the non-Gaussianity of the 6Ti, field is less than the mass density field of 
hydrogen clouds. This is because high-density areas correspond to lower ffji regions. The 
non-Gaussianity is weakened when transferring the mass field to the ffji field by equation 
(11). Moreover, equation (12) leads to further weakening of the non-Gaussianity given 
by high-density peaks. Therefore, although its PDFs (Fig. 7) are long tailed, the non- 
Gaussianity of the 6Tb field is not as strong as that of an LN field. 

High-order moments are more sensitive to the tails of 5T5 PDFs, and therefore, the 
statistics {{6Tb - {6Tb)f'')'^''^/{{6Tb - {6Tb)f) and {{6T^jy/''/{{6Tb)) are less dependent on 
the lower limit of observable 5T5. If we use only the observable events with 6Th > 1 mK, 
the high-order moments given by the truncated PDFs [Fig. 7] will not be changed very 
much. Moreover, the statistics {{6Tb-{6Tb)f'^y/'^/{{6Tb-{6Tb))^) and {{6T^^^)^I'' / {{6Tb)) 
are defined by ratios of 6Tb. Hence, they are less dependent on mean fui- The statistical 
behavior shown in Figure 11 is weakly dependent on ionizing photon parameter J in equation 
(11). 



4.5. Scale-scale correlations of the WFCs 

The so-called scale-scale correlations of a random field measures the correlation between 
the fiuctuations on different scales (Pando et al. 1998). For 6Tb field, this statistics is defined 
as 

where p is an even integer. The notation [...] in equation (26) denotes the integer part of 
the quantity. Because the spatial range of the cell (j, [1/2]) is the same as that of the two 
cells (j + 1, /) and (j + 1, / + 1), Cj'^ measures the correlation between the fiuctuations on 
scales j and j -|- 1 at the same physical area. For Gaussian noise, the variables ej^i on scales 

i and j + 1 are uncorrected, and hence we have (e?rw2i^+i i) = (^fi/2i)(^+i 1)' '^^ ^fj+i ~ ^' 
Therefore, the scale-scale correlation is effective for drawing a non-Gaussian signal from a 
Gaussian background, even when the variance of the noise is comparable to the signal (Feng 
et al. 2000; Feng & Fang 2000). 



^p,p _ \''j,[l/2fj+l,l/ , 

^j - Tip rr=P T' l^*^^ 



Figure 12 presents the C^'^ versus j of the 6Tb field with parameters (Tq, rj, z) = (1.3 x 
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W'^K, 2, 7) and (0.8 x WK, 2, 10) and p = 2, 4 and 8. Figure 12 shows once again that 
the statistical resuhs are less dependent on the temperature parameters considered. For all 
cases of p = 2, 4 and 8, the ST^ fields are substantially scale-scale correlated on scales j > 3, 
corresponding to length scale 6.25 h~^ Mpc and angular scale 3'. 52 at 2; = 7 or 3'. 28 at 
z = 10. At scale j = 8, or 0.20 /i^^Mpc, Cj'^ is as large as 10^, which would be very useful 
for distinguishing the signal of 6Tb{zo) from noise. 

From Figure 12 we see that the scale-scale correlation is remarkably increasing from 
p = 2 to p = 4, but not so from p = 4 to p = 8. This is again because the PDFs of 6Ti, 
are long tailed, but not very long. This can strongly affect statistics of order n < 5, but 
not n > 5. We calculate also the scale-scale correlation for noise-contaminated samples. It 
indeed shows that the Gaussian noise can be filtered out with the scale-scale correlations. 



5. Discussion and Conclusion 

We have studied the brightness temperature excess STb of the redshifted 21 cm emission 
of hydrogen clouds in the epochz^ > z > zgp- According to the LN model, the mean 
temperature of hydrogen clouds does not undergo strong evolution and is generally Tq ~ 10^ 
K. Therefore, the 21 cm emission from this epoch possesses some statistical features weakly 
dependent on the details of this epoch. They are as follows: 

1. The random field of 6Tb is substantially non-Gaussian. It consists of spikes with high 
STf, and a low STb area between the spikes. 

2. The mean of 5Tb is ~ 2 mK, and the variance ((5Tf,)^)^/^ is of the order of 10 mK on 
about arcminute scales, while the spikes can be as high as 40 niK. 

3. The one-point distributions of either 6Tb or 6Tb^(j^i) have long tails approximately 
following power laws. 

4. The nth-order moment of STb is quickly increasing with n when n < 5, and slowly 
when n > 5. It is very different from a Gaussian field and an LN field. 

5. The scale-scale correlation is significant for all scales, which is effective for washing 
out Gaussian noise contamination. 

All these results show that although the clustering of hydrogen gas is suppressed at 
z < Zr, the 21 cm emission fields are significantly non-Gaussian since Zr- 

To obtain the above results, we used the fraction of neutral hydrogen HI given by the 
LN model (Liu et al. 2004). Actually, the evolution of the HI fraction in Zr > z > zqp is 
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not strongly constrained by the current observations. The optical depth Te is sensitive to Zr 
and the mean of the HI fraction in Zj. > z > zgp, but not to the details of the evolution 
of the HI fraction. Different evolution models of the HI and HII fractions can reasonably 
explain the observed Tg. However, among the above results, the last two points (high-order 
moments and scale-scale correlations) are based on statistics given by the ratios of STh. They 
are probably less dependent on the HI or HII fractions. 

We show that the non-Gaussian properties of the STh field comes from the non-Gaussianity 
of the mass field at high redshifts. For a given redshift z, the STf, field of the redshifted 21 
cm emission is dependent only on T^ at that redshift, but less sensitive to the details of 
the redshift evolution of Tm- On the other hand, in the epoch Zr > z > zgp, Tm is in the 
range 10^ — 10^ K. Therefore, the non-Gaussian features revealed in this paper may help 
in, or even be necessary for, unambiguously identifying the 21 cm emission from the epoch 
Zr > z > Zgp- 
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Fig. 1. — Comoving Jeans length xj as a function of the cosmic scale factor 1/(1 + z). 
In the period Zr < z < zgp, the mean temperature of hydrogen is assumed to be T^ = 
To[(l + z)/{l + Zr)]"^, with To = 1.3 x 10^ K (solid line) and 0.8x10^ K (dotted line). 

Fig. 2. — Variance o"o of linear perturbations as a function of 1/(1 + z), which is calculated 
with the xj given by Fig. 1. 

Fig. 3. — Temperature evolution in the three epochs: (1) z > Zr = 18, T = 4[(1 + z)/16]^ 
K; (2) Zr> z> zgp = 7,T = To[(l + z)/{l + Zr)f, with Tq = 1.3 x 10^ K (solid line), and 
0.8x10"^ K (dotted line); (3) z < zgp, T ~ Tq. The increase oi T a.t z = 3.3 is due to He 
reionization. 

Fig. 4. — Realization of hydrogen density fields at redshifts 7, 10 and 15. The densities are 
in units of picM- 

Fig. 5. — Realization of the fields of Tf,, T21 and STf,. The relevant parameters {To,r],z) are 
taken to be (1.3 x 10^ K, 2, 7). 

Fig. 6. — Cumulative {5T^)^'^{> p) mK for field with density > p as a function of p. The 
parameters {TQ,r],z) are taken to be (1.3 x 10*^ K, 2, 7). 

Fig. 7. — One-point distributions of 57b^(j,«) for (To,ri,z) = (1.3 x 10'^ K, 2, 7) (left panels) 
and (0.8 x 10^ K, 2, 10) (right panels). The smoothing scales are j = 6, 8, and 10. 

Fig. 8. — One-point distributions of STh^i^j^i^ for {TQ,r],z) = (1.3 x 10'^ K, 2, 7) (left panels) 
and (0.8 x 10^ K, 2, 10) (right panels). The smoothing scales are j = 6, 8, and 10. 

Fig. 9. — Two-mode correlation function ln(5Tfe_Qy)(5Tb (j;/)) vs. Inr {h^^ Mpc) for parameters 
(To,ri,z) = (1.3 X 10'^ K, 2, 7). The smoothing scales are j = 8, 10, and 12. 

Fig. 10. — DWT power spectrum Pj of the STb field for parameters (Tq, rj, z) = (1.3 x 10^ K, 
2, 7) and (0.8 x 10^ K, 2, 10). 

Fig. 11.— The n-dependence of (a) {{6Thj - {6Thj)f''y/''/{{6Tbj - {SThj)f) and (b) 
{{ST^jy/'' / {{5Tb)) for parameters {To,ri,z) = (1.3 x 10^ K, 2, 7) (left panels) and (0.8 x 10^ 
K, 2, 10) (right panels) 

Fig. 12. — Scale-scale correlation Cj'^ vs. j for parameters {To,ri,z) = (1.3 x 10^ K, 2, 7) 
(left panels) and (0.8 x 10*^ K, 2, 10) (right panels); p is taken to be 2, 4, and 8. 
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